% Function used to obtain moments without stacking each variable into one
% column with all (nf x T) observations.
function M = getmoments(DATAJ,nf,nx,T)

% 1) Corr(I/K(t),I/K(t-1))
% 2) Corr(I/K(t),H/L(t-1))
% 3) Corr(I/K(t),Y_growth(t-1))
% 4) Corr(H/L(t),I/K(t-1)) 
% 5) Corr(H/L(t),H/L(t-1))
% 6) Corr(H/L(t),Y_growth(t-1))
% 7) Corr(Y_growth(t),I/K(t-1))
% 8) Corr(Y_growth(t),H/L(t-1))
% 9) Corr(Y_growth(t),Y_growth(t-1))
% 10) Corr(I/K(t),I/K(t-2))
% 11) Corr(I/K(t),H/L(t-2))
% 12) Corr(I/K(t),Y_growth(t-2))
% 13) Corr(H/L(t),I/K(t-2))
% 14) Corr(H/L(t),H/L(t-2))
% 15) Corr(H/L(t),Y_growth(t-2))
% 16) Corr(Y_growth(t),I/K(t-2))
% 17) Corr(Y_growth(t),H/L(t-2))
% 18) Corr(Y_growth(t),Y_growth(t-2))
% 19) Corr(I/K(t),H/L(t))
% 20) Corr(I/K(t),Y_growth(t))
% 21) Corr(H/L(t),Y_growth(t))
% 22) Mean(I/K)
% 23) Std(I/K)
% 24) Mean(H/L)
% 25) Std(H/L)
% 26) Std(Y_growth)
% 27) Sample proportion: Negative investments
% 28) Sample proportion: Zero investments
% 29) Sample proportion: Positive investments
% 30) Sample proportion: Negative gross hirings
% 31) Sample proportion: Zero gross hirings
% 32) Sample proportion: Positive gross hirings
% 33) Sample proportion: Negative net hirings
% 34) Sample proportion: Zero net hirings
% 35) Sample proportion: Positive net hirings
% 36) Sample proportion: Negative investments & negative gross hirings
% 37) Sample proportion: Negative investments & zero gross hirings
% 38) Sample proportion: Negative investments & positive gross hirings
% 39) Sample proportion: Zero investments & negative gross hirings
% 40) Sample proportion: Zero investments & zero gross hirings
% 41) Sample proportion: Zero investments & positive gross hirings
% 42) Sample proportion: Positive investments & negative gross hirings
% 43) Sample proportion: Positive investments & zero gross hirings
% 44) Sample proportion: Positive investments & positive gross hirings
% 45) Sample proportion: Negative investments & negative net hirings
% 46) Sample proportion: Negative investments & zero net hirings
% 47) Sample proportion: Negative investments & positive net hirings
% 48) Sample proportion: Zero investments & negative net hirings
% 49) Sample proportion: Zero investments & zero net hirings
% 50) Sample proportion: Zero investments & positive net hirings
% 51) Sample proportion: Positive investments & negative net hirings
% 52) Sample proportion: Positive investments & zero net hirings
% 53) Sample proportion: Positive investments & positive net hirings
% 54) Sample proportion: Non-zero investments & non-zero gross hirings
% 55) Sample proportion: Non-zero investments & zero gross hirings
% 56) Sample proportion: Zero investments & non-zero gross hirings
% 57) Sample proportion: Zero investments & zero gross hirings
% 58) Sample proportion: Non-zero investments & non-zero net hirings
% 59) Sample proportion: Non-zero investments & zero net hirings
% 60) Sample proportion: Zero investments & non-zero net hirings
% 61) Sample proportion: Zero investments & zero net hirings
% 62) First order serial correlation, capital.
% 63) Second order serial correlation, capital.
% 63) First order serial correlation, labor.
% 65) Second order serial correlation, labor.
% 66) Sample cross correlation, log(capital) and log(sales)
% 67) Sample cross correlation, log(labor) and log(sales)
% 68) Sample cross correlatoin, log(capital) and log(sales)
% 69) Sample skewness in distribution of investment rates.
% 70) Corr(I/K(t), log(Y/K(t)))
% 71) Corr(I/K(t), log(Y/K(t-1)))
% 72) Corr(I/K(t), log(Y/K(t-2)))
% 73) Corr(H/L(t), log(Y/L(t)))
% 74) Corr(H/L(t), log(Y/L(t-1)))
% 75) Corr(H/L(t), log(Y/L(t-2)))
% 76) Corr(H/L(t), log(Y/K(t)))
% 77) Corr(H/L(t), log(Y/K(t-1)))
% 78) Corr(H/L(t), log(Y/K(t-2)))
% 79) Corr(Ygr(t), log(Y/K(t)))
% 80) Corr(Ygr(t), log(Y/K(t-1)))
% 81) Corr(Ygr(t), log(Y/K(t-2)))
% 82) Corr(Ygr(t), log(Y/L(t)))
% 83) Corr(Ygr(t), log(Y/L(t-1)))
% 84) Corr(Ygr(t), log(Y/L(t-2)))
% 85) Corr(log(Y/K(t)), log(Y/K(t-1)))
% 86) Corr(log(Y/K(t)), log(Y/K(t-2)))
% 87) Corr(log(Y/L(t)), log(Y/L(t-1)))
% 88) Corr(log(Y/L(t)), log(Y/L(t-2)))
% 89) Sample skewness in distribution of hiring rates.
% 90) Sample proportion of negative investments (i<-0.02)
% 91) Sample proportion of zero investments     (-0.02<i<0.02)
% 92) Sample proportion of zero gross hirings   (-0.02<h<0.02)
% 93) Sample proportion of zero net hirings     (-0.02<dL<0.02)
% 94) Sample proportion of simultaneous zeros, gross hirings (-0.02<i<0.02 & -0.02<h<0.02)
% 95) Sample proportion of simultaneous zeros, net hirings   (-0.02<i<0.02 & -0.02<dL<0.02)   

% Determine some inaction definitions:
IKinaction = [-0.03 0.03];
HLgrossinaction = [-0.03 0.03];
HLnetinaction = [0.07 0.13];

% Declare empty moment vector:
M = [];

% Define matrix of investment rates, dimension (T x nf):   
IK = zeros(T,nf);
for i = 1 : nf;
  IK(:,i) = DATAJ(:,nx*(i-1)+6);
end

% Define matrix of hiring rates, dimension (T x nf):   
HL = zeros(T,nf);
for i = 1 : nf;
  HL(:,i) = DATAJ(:,nx*(i-1)+7);
end

% Define matrix of capital input, dimension (T x nf):
K = zeros(T,nf);
for i = 1 : nf;
  K(:,i) = DATAJ(:,nx*(i-1)+4);
end

L = zeros(T,nf);
% Define matrix of labor input, dimension (T x nf):
for i = 1 : nf;
  L(:,i) = DATAJ(:,nx*(i-1)+5);
end

% Define matrix of sales values, dimension (T x nf):  
Y = zeros(T,nf);
for i = 1 : nf;
  Y(:,i) = DATAJ(:,nx*(i-1)+3);
end

% Define matrix of sales growth rates, dimension (T-1 x nf):
Ygr=zeros(T-1,nf);
for i = 1 : nf;
  Ygr(:,i) = DATAJ(2:end,nx*(i-1)+3)./DATAJ(1:end-1,nx*(i-1)+3)-1;
end

YK = Y./K;
yk = log(YK);
YL = Y./L;
yl = log(YL);

% Fill in some moments:
M(1,22) = mean(reshape(IK,T*nf,1));
M(1,23) = std(reshape(IK,T*nf,1));
M(1,24) = mean(reshape(HL,T*nf,1));
M(1,25) = std(reshape(HL,T*nf,1));
M(1,26) = std(reshape(Ygr,(T-1)*nf,1));
M(1,69) = skewness(reshape(IK,T*nf,1));
M(1,89) = skewness(reshape(HL,T*nf,1));
M(1,27) = mean(reshape(IK<IKinaction(1),T*nf,1));
M(1,28) = mean(reshape((IK>=IKinaction(1)).*(IK<=IKinaction(2)),T*nf,1));
M(1,29) = mean(reshape(IK>IKinaction(2),T*nf,1));
M(1,30) = mean(reshape(HL<HLgrossinaction(1),T*nf,1));
M(1,31) = mean(reshape((HL>=HLgrossinaction(1)).*(HL<=HLgrossinaction(2)),T*nf,1));
M(1,32) = mean(reshape(HL>HLgrossinaction(2),T*nf,1));
M(1,33) = mean(reshape(HL<HLnetinaction(1),T*nf,1));
M(1,34) = mean(reshape((HL>=HLnetinaction(1)).*(HL<=HLnetinaction(2)),T*nf,1));
M(1,35) = mean(reshape(HL>HLnetinaction(2),T*nf,1));
M(1,36) = mean(reshape((IK<IKinaction(1)).*(HL<HLgrossinaction(1)),T*nf,1));
M(1,37) = mean(reshape((IK<IKinaction(1)).*(HL>=HLgrossinaction(1)).*(HL<=HLgrossinaction(2)),T*nf,1));
M(1,38) = mean(reshape((IK<IKinaction(1)).*(HL>HLgrossinaction(2)),T*nf,1));
M(1,39) = mean(reshape((IK>=IKinaction(1)).*(IK<=IKinaction(2)).*(HL<HLgrossinaction(1)),T*nf,1));
M(1,40) = mean(reshape((IK>=IKinaction(1)).*(IK<=IKinaction(2)).*(HL>=HLgrossinaction(1)).*(HL<=HLgrossinaction(2)),T*nf,1));
M(1,41) = mean(reshape((IK>=IKinaction(1)).*(IK<=IKinaction(2)).*(HL>HLgrossinaction(2)),T*nf,1));
M(1,42) = mean(reshape((IK>IKinaction(2)).*(HL<HLgrossinaction(1)),T*nf,1));
M(1,43) = mean(reshape((IK>IKinaction(2)).*(HL>=HLgrossinaction(1)).*(HL<=HLgrossinaction(2)),T*nf,1));
M(1,44) = mean(reshape((IK>IKinaction(2)).*(HL>HLgrossinaction(2)),T*nf,1));
M(1,45) = mean(reshape((IK<IKinaction(1)).*(HL<HLnetinaction(1)),T*nf,1));
M(1,46) = mean(reshape((IK<IKinaction(1)).*(HL>=HLnetinaction(1)).*(HL<=HLnetinaction(2)),T*nf,1));
M(1,47) = mean(reshape((IK<IKinaction(1)).*(HL>HLnetinaction(2)),T*nf,1));
M(1,48) = mean(reshape((IK>=IKinaction(1)).*(IK<=IKinaction(2)).*(HL<HLnetinaction(1)),T*nf,1));
M(1,49) = mean(reshape((IK>=IKinaction(1)).*(IK<=IKinaction(2)).*(HL>=HLnetinaction(1)).*(HL<=HLnetinaction(2)),T*nf,1));
M(1,50) = mean(reshape((IK>=IKinaction(1)).*(IK<=IKinaction(2)).*(HL>HLnetinaction(2)),T*nf,1));
M(1,51) = mean(reshape((IK>IKinaction(2)).*(HL<HLnetinaction(1)),T*nf,1));
M(1,52) = mean(reshape((IK>IKinaction(2)).*(HL>=HLnetinaction(1)).*(HL<=HLnetinaction(2)),T*nf,1));
M(1,53) = mean(reshape((IK>IKinaction(2)).*(HL>HLnetinaction(2)),T*nf,1));
M(1,54) = mean(reshape(((IK<IKinaction(1))+(IK>IKinaction(2))).*((HL<HLgrossinaction(1))+(HL>HLgrossinaction(2))),T*nf,1));
M(1,55) = mean(reshape(((IK<IKinaction(1))+(IK>IKinaction(2))).*(HL>=HLgrossinaction(1)).*(HL<=HLgrossinaction(2)),T*nf,1));
M(1,56) = mean(reshape((IK>=IKinaction(1)).*(IK<=IKinaction(2)).*((HL<HLgrossinaction(1))+(HL>HLgrossinaction(2))),T*nf,1));
M(1,57) = mean(reshape((IK>=IKinaction(1)).*(IK<=IKinaction(2)).*(HL>=HLgrossinaction(1)).*(HL<=HLgrossinaction(2)),T*nf,1));
M(1,58) = mean(reshape(((IK<IKinaction(1))+(IK>IKinaction(2))).*((HL<HLnetinaction(2))+HL>HLnetinaction(2)),T*nf,1));
M(1,59) = mean(reshape(((IK<IKinaction(1))+(IK>IKinaction(2))).*(HL>=HLnetinaction(1)).*(HL<=HLnetinaction(2)),T*nf,1));
M(1,60) = mean(reshape((IK>=IKinaction(1)).*(IK<=IKinaction(2)).*((HL<HLnetinaction(1))+(HL>HLnetinaction(2))),T*nf,1));
M(1,61) = mean(reshape((IK>=IKinaction(1)).*(IK<=IKinaction(2)).*(HL>=HLnetinaction(1)).*(HL<=HLnetinaction(2)),T*nf,1));

% Determine some inaction definitions:
IKinaction = [-0.02 0.02];
HLgrossinaction = [-0.02 0.02];
HLnetinaction = [0.08 0.12];

M(1,90) = mean(reshape(IK<IKinaction(1),T*nf,1));
M(1,91) = mean(reshape((IK>=IKinaction(1)).*(IK<=IKinaction(2)),T*nf,1));
M(1,92) = mean(reshape((HL>=HLgrossinaction(1)).*(HL<=HLgrossinaction(2)),T*nf,1));
M(1,93) = mean(reshape((HL>=HLnetinaction(1)).*(HL<=HLnetinaction(2)),T*nf,1));
M(1,94) = mean(reshape((IK>=IKinaction(1)).*(IK<=IKinaction(2)).*(HL>=HLgrossinaction(1)).*(HL<=HLgrossinaction(2)),T*nf,1));
M(1,95) = mean(reshape((IK>=IKinaction(1)).*(IK<=IKinaction(2)).*(HL>=HLnetinaction(1)).*(HL<=HLnetinaction(2)),T*nf,1));

% Get coefficents of serial and cross correlation:
C = autocorr(reshape(IK,size(IK,1)*size(IK,2),1),2,[]);
M(1,1)  = C(2);
M(1,10) = C(3);
C = autocorr(reshape(HL,size(HL,1)*size(HL,2),1),2,[]);
M(1,5)  = C(2);
M(1,14) = C(3);
C = autocorr(reshape(Ygr,size(Ygr,1)*size(Ygr,2),1),2,[]);
M(1,9)  = C(2);
M(1,18) = C(3);
C = autocorr(reshape(K,size(K,1)*size(K,2),1),2,[]);
M(1,62) = C(2);
M(1,63) = C(3);
C = autocorr(reshape(L,size(L,1)*size(L,2),1),2,[]);
M(1,64) = C(2);
M(1,65) = C(3);
C = crosscorr(reshape(IK,size(IK,1)*size(IK,2),1),reshape(HL,size(HL,1)*size(HL,2),1),2,[]);
M(1,2)  = C(2);
M(1,11) = C(1);
M(1,4)  = C(4);
M(1,13) = C(5);
M(1,19) = C(3);
C = crosscorr(reshape(IK(2:end,:),(size(IK,1)-1)*size(IK,2),1),reshape(Ygr,size(Ygr,1)*size(Ygr,2),1),2,[]);
M(1,3)  = C(2);
M(1,12) = C(1);
M(1,7)  = C(4);
M(1,16) = C(5);
M(1,20) = C(3);
C = crosscorr(reshape(HL(2:end,:),(size(HL,1)-1)*size(HL,2),1),reshape(Ygr,size(Ygr,1)*size(Ygr,2),1),2,[]);
M(1,6)  = C(2);
M(1,15) = C(1);
M(1,8)  = C(4);
M(1,17) = C(5);
M(1,21) = C(3);
C = crosscorr(reshape(log(K),size(K,1)*size(K,2),1),reshape(log(Y),size(Y,1)*size(Y,2),1),2,[]);
M(1,66) = C(3);
C = crosscorr(reshape(log(L),size(L,1)*size(L,2),1),reshape(log(Y),size(Y,1)*size(Y,2),1),2,[]);
M(1,67) = C(3);
C = crosscorr(reshape(log(K),size(K,1)*size(K,2),1),reshape(log(L),size(L,1)*size(L,2),1),2,[]);
M(1,68) = C(3);
C = crosscorr(reshape(IK,size(IK,1)*size(IK,2),1),reshape(yk,size(yk,1)*size(yk,2),1),2,[]);
M(1,70) = C(3);
M(1,71) = C(2);
M(1,72) = C(1);
C = crosscorr(reshape(HL,size(HL,1)*size(HL,2),1),reshape(yl,size(yl,1)*size(yl,2),1),2,[]);
M(1,73) = C(3);
M(1,74) = C(2);
M(1,75) = C(1);
C = crosscorr(reshape(HL,size(HL,1)*size(HL,2),1),reshape(yk,size(yk,1)*size(yk,2),1),2,[]);
M(1,76) = C(3);
M(1,77) = C(2);
M(1,78) = C(1);
C = crosscorr(reshape(Ygr,size(Ygr,1)*size(Ygr,2),1),reshape(yk(2:end,:),(size(yk,1)-1)*size(yk,2),1),2,[]);
M(1,79) = C(3);
M(1,80) = C(2); 
M(1,81) = C(1);
C = crosscorr(reshape(Ygr,size(Ygr,1)*size(Ygr,2),1),reshape(yl(2:end,:),(size(yl,1)-1)*size(yl,2),1),2,[]);
M(1,82) = C(3);
M(1,83) = C(2);
M(1,84) = C(1);
C = autocorr(reshape(yk,size(yk,1)*size(yk,2),1),2,[]);
M(1,85) = C(2);
M(1,86) = C(3); 
C = autocorr(reshape(yl,size(yl,1)*size(yl,2),1),2,[]);
M(1,87) = C(2);
M(1,88) = C(3);

